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I. INTRODUCTION 


Our laboratory is presently investigating possible instrumentation for geo- 
chemical exploration of lunar or planetary surfaces for both manned and unmanned 
missions. The studies presently going forward will hopefully be used to define a 
total system having as its objectives: (1) obtaining operational numbers to help 
an astronaut in sample selection, and (2) obtaining geochemical information for 
mapping and preliminary surface analysis for manned and unmanned missions. 

During the NASA 1965 Summer Conference on Lunar Exploration and 
Science, a considerable amount of attention was given by the Geochemistry Work- 
ing Group to problems of surface exploration. While it was the consensus that 
the lunar environment was not particularly conducive to performing detailed geo- 
chemical analysis, there was, however, a great need for diagnostic tools to help 
the astronaut in selecting samples to be returned to earthside laboratories for 
comprehensive studies. 

Consideration of diagnostic tools in relationship to the mission’s objectives 
leads to a number of requirements that must be satisfied. These requirements 
range from geochemical parameters to be tested, optimum instruments to be 
employed and the vital problem of data acquisition and processing in very real 
time. Additionally, although the proposed instrumentation would primarily enable 
the astronaut to make decisions about sample selection, the design is sufficiently 
flexible to be useful for either manned or unmanned missions. 

There are certain obvious characteristics that such instrumentation must 
meet such as small size and weight, reliability, stability and specificity. Be- 
cause some of these requirements are at least in part conflicting, one must- 
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accept some tradeoffs. In the paper that follows one such instrument will be 
described as well as a computer program for data reduction. 

A. Instrumentation 

Obviously, instrumental and system design is strictly determined by the 
mission profile. If, for example, one established a criterion that the probe is 
needed to help an astronaut in sample selection, by permitting him to perform 
relatively simple "in situ" analysis, then "non-dispersive" x-ray emission* 
analysis becomes attractive. Such a device lends itself to portability as well as 
simple and rugged construction having minimal space and power requirements. 

In fact, such instrumentation has been proposed by a number of investigators 
such as Metzger, et al. , Sellers and Ziegler , and Kartunnen, et al. , Trombka 
and Adler 4 . One of the most desirable features of such instrumentation is that 
x-ray spectra can be effectively excited by means of radioactive isotopic sources. 
The use of this type excitation has been discussed in a number of publications by 
Rhodes 5 , Cameron and Rhodes 6 , Robert 7 , Friedman 8 , Robert and Martinelli 9 , 
Inamura, et al. 10 . Most applications have involved the use of gamma, beta or 
bremsstrahlung sources but only recently has a systematic study of alpha exci- 
tation been undertaken.** 

The advantages of alpha particle excitation will be discussed in additional 
detail below. 

In the program to be described attention has been directed to the common 
rock forming elements such as Mg, Al, Si, K, Ca and Fe. The characteristic K 


*The term non-dispersive is actually a misnomer. This expression which is in common usage 
means simply pulse height analysis as opposed to Bragg-Crystal dispersion. 

**See References 2, 7, 8, 10. 
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x-rays emitted by the majority of these elements may be classed as soft x-rays 
and, with the exception of Fe, range from 3A to approximately 12A. In this range 
one encounters problems of excessive absorption in the detector windows as 
well as low fluorescence yields. The instrument therefore requires an efficient 
design with emphasis on optimum excitation and geometry. It is also essential 
to use the thinnest detector windows compatible with mechanical strength and low 
counter gas leakage. An additional but very important requirement is to obtain 
the best possible detector resolution. 

Figure 1 shows the rather simple instrumentation used in the laboratory* 
for analysis. The apparatus consists of a radioactive source, a means for hold- 
ing the sample, a proportional detector and an inlet for introducing helium into 
the chamber. Although a modest vacuum is desirable, helium is used at present 
in order to preserve the integrity of the alpha source window. 

The samples to be analyzed may be introduced as powders, briquets or 
small flat sections of rock. The various components of the instrumentation will 
now be examined more closely. 

B. Excitation Sources 

The obvious advantage in using radio-isotopes for exciting x-rays is that the 
source is small, requires no power and is the ultimate in reliability. It has al- 
ready been amply demonstrated that radioactive sources of relatively low output 
will generate very adequate x-ray spectra in time periods as short as minutes. 

Various types of isotopes have been found useful for x-ray emission analyses. 
These are beta sources, bremsstrahlung sources, K capture x-ray sources, and 

*This is not proposed as a flight instrument but is rather a breadboard for laboratory evaluation of 
the technique. 
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alpha emitters. It has been shown by Sellers and Ziegler that alpha sources 
are particularly useful for light element analysis. These advantages can be 
tabulated as follows: 

1. The cross section for ionization of the K x-rays increases approximately 
as l/Z 12 . Thus the rapidly increasing ionization cross section with de- 
creasing atomic number more than compensates for the decreasing 
fluorescence yield. 

2. In comparison to electron excitation the alpha bremsstrahlung produced 
continuum is reduced by a factor of (m/M) where m and M are, re- 
spectively, the masses of the electron and alpha particles. As a conse- 

6 7 

quence the continuum is reduced by a factor of about 10 ot 10 . 

3. Although alpha emission is usually associated with gamma ray emission 
it is possible to select nuclides where the gamma emission is both a 
small fraction of the alpha emission and the gamma rays are sufficiently 
removed from the excited x-rays to be easily discriminated against. 

4. It is possible to obtain sources of high specific activity emitting very 

210 

energetic alpha particles having energies of 4 to 6 mev such as Po 
and Cm 242 . 

While the advantages listed above make alpha emitters attractive sources 
for light element analysis there are a number of problems encountered in 
practice. Alpha sources are a serious potential health hazard and must be 
handled with extreme care. Because of the phenomenon of aggregate recoil, one 
should work either with sealed sources or in adequately ventilated hoods or 
glove boxes. As yet it is difficult to obtain sealed sources that can be operated 
in vacuum without using windows of such thickness as to seriously degrade the 
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energy of the emitted alpha particles. In the work to be described here measure- 
ments were made in helium at atmospheric pressure. The source holder is shown 

242 

in Figure 2. The Cm source was approximately 15 me deposited by coating 
on a stainless steel foil. The window in this instance was 0.0005 inch aluminum. 
Although no direct measurement of the energy distribution of the alpha particles 
was made, it was estimated that the effect of the window and helium atmosphere 
reduced the 6 mev alphas about 3 mev with an associated broadening of the 
energy peak. 

Although it is true that alpha bombardment excites virtually no bremsstrah- 

lung, one must expect some background. This background may come from the 

source container as well as x-rays excited by processes such as internally 

converted gamma rays and due to contaminating radioactive species. Figure 3 

242 

shows the background radiation from the Cm source scattered from a plastic 
sample. One can identify at the very least an AIK line from the source window, 
an FeK line from the substrate on which the Cm is deposited, a W line from 
the hevi-met (sintered W) container, as well as lines in the 13-18 kev region 
which may be attributed to the daughter product Pu as well as other contami- 
nants. 

C. Detector 

A major problem in light element analysis is associated with detection. The 
ideal would be a windowless detector with high efficiency and resolution. This is 
unfortunately beyond the state of the art at this time, particularly in the region 
of from 7 to 1 kev. The best detector available at this time is the thin window 
proportional counter, either flow or sealed. In the case of flow counters one can 
work easily with 0.001 inch Be or 0.00025 inch mylar. Figure 4 shows a typical 
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efficiency curve for a 2 cm detector using an Argon-methane mixture. This 
curve represents the product of window transmission and gas absorption. It can 
be seen that counter efficiency drops to about 10% at 11 A. When one combines 
this with phenomenon with the low fluorescence yield for the low Z elements, 
then the high ionization cross section of the alpha particles becomes particularly 
significant. The detector resolution used in this study was of the order of 17% 
for Fe 55 . 


D. Electronics 

This electronic arrangement was conventional consisting of a charge sensi- 
tive preamplifier feeding a high quality multi-channel analyzer. Although 
greater capabilities were available, 256 channels of the analyzer were found to 
be sufficient. Data readout was on perforated tapes which were then converted 
to punch cards for computer processing. 

D. DATA ANALYSIS 
A. Analytical Approach 

In the previous section, the design of the X-ray emission system was shown 
to be greatly influenced by the mission objectives. Similarly, the analytical 
method is also closely related to these objectives. Methods can be developed 
for analysis of the measured spectra which at one extreme can compare spectral 
patterns and at the other extreme decompose the observed spectrum into mono- 
elemental components for detailed analysis. The simpler pattern recognition 
techniques are now being investigated for their possible application to problems 
involving sample selection, but since in almost all cases the more detailed 
analysis will be required, we will address ourselves to this problem. 
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The problem of determining qualitative and quantitative information from a 
non-dispersive measurement of the x-ray emission can be divided into two parts. 
First, the differential x-ray energy spectrum (true energy spectrum incident 
upon the proportional counter) must be determined from the measured pulse 
height spectrum. The second step is then to deduce the chemical composition 
of the irradiated sample from the differential energy x-ray spectrum. The first 
step will be described in detail in this paper, and possible methods for perform- 
ing the second step will be discussed by example at the end of this section. 

B. Nature of the Pulse Height Spectrum 
We shall begin by defining what is meant by a pulse height spectrum. Con- 
sider the case of an ideal spectrometer looking at the Fe K spectrum, consist- 
ing of a number of monochromatic energies. We would observe a spectral 
distribution consisting of a number of discrete lines as shown in Figure 5. If, in 
practice, the x-ray energies are examined by means of a proportional counter 
detector and a pulse height analyzer, the spectrum becomes a continuous distri- 
bution as shown in Figure 5. Here we see the effect of instrumental smearing 
due to such factors as resolution, a statistical fluctuation in the detection 
process, detector geometry, electronic noise, etc. In many cases the spectrum 
is further complicated by an escape peak phenomenon due to partial absorption 
of the incident energies. It is obvious that any program of reducing a pulse 
height distribution to pulse height energies must concern itself with and remove 
these instrumental factors. This necessary approach will now be discussed in 
some detail. 
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It has been observed that the x-ray photopeaks show a nearly Gaussian 
distribution and that the resolution of the peaks is a function of energy, the 
resolution R being defined as 


where W 1/2 is the width of the photopeak at half the maximum amplitude, and 
P m is the position of the pulse height maximum, where the width and peak are 
expressed in the same units (energy or channel window). The resolution is in- 
versly proportional to some power of the pulse height or energy 


R ~ E -n 


( 2 ) 


Theoretically 11 n = 0.5 although experimentally n is usually greater or equal 
to 0.5. In the cases we have examined, n has varied from 0.5 to 0.67 and ap- 
pears to be affected by the resolution of a given detector. 

If one irradiates a mixture of elements with alpha particles and records the 
pulse height spectrum, the total observed poly elemental spectrum will be made 
up of a summation of the monoelemental components as shown in Figure 6. 

If /Oj is taken as the measured count in channel i, then 

* = L <3 > 

where Cix. is the total number of counts in channel i due to the A.th element. 
The pulse height spectra of the monoelemental components can be normalized 
as follows: 



Pk A i> 


( 4 ) 
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Pulse height spectrum of a basalt 
measured with a proportional counter 
filled with 90% A and 10% CH . The 
following can be seen: 

1 . The monoelemental components 

2. The envelope due to the sum of 
the monoelemental components 

3. The measured pulse height spectrur 
with the statistical variations 
included 





where is the normalized number of counts in channel i due to the element 
k, and /3^ is the relative intensity of the k component, then 



Another component of the measured pulse height spectra which must be 
considered is the background. In the case of alpha induced fluorescence there 
are two separate sources of radiation contributing to the background spectrum 
(see Figures 3 and 6). As indicated above, one component can be attributed to 
the natural background and to gamma rays from the source striking the detector. 
The second component can be attributed to the x-rays emanating from the source 
being scattered by the sample. We have found that in this instance the major 
contribution to the background seems to be due to this second scattered com- 
ponent. The amount of scatter will depend on the average atomic number of the 
sample and thus changes from sample to sample. We have found further from 
our measurements that while the shape of the background pulse height spectrum 
does not change significantly from sample to sample the intensity does. Further 
study of the shape of the background as a function of composition is underway. 
The problem of subtracting background is not a simple one. Our approach to 
this problem will be described in detail, in the final section of this paper where' 
actual spectral analyses will be presented. Figure 6 shows the measured pulse 
height spectrum of a basalt sample with the background included. Equation (3a) 
should include the background component B., where B i (see Figure 3) is the 
intensity of the background in channel i : 

♦ B r < 3b > 

\ * 
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C. Formulation of the Linear Least Square Analysis 
1. Matrix Equations: If equation (3b) were truly an equality the relative 
intensities /3^ could easily be determined by a simple matrix inversion of (3b) 
that is 


yS = (A -*)p, (4) 

where f3 is a vector of the relative intensities, 

A is a m + n matrix of the monoelemental components, and 
p is a vector of the measured polyelemental spectrum minus background 
(P'i - Bj). 

The problem in performing this inversion arises when one considers the 
statistical variance in the measurement Aj^ , and p .. Because of this, the 
equality (3b) is not true and a unique solution to (4) does not exist. Equation (4) 
is overdetermined, therefore an infinite set of equally good solutions is possible 
and there are no criteria for selecting the best one. As a consequence of this 
and because of the statistical variance involved only the most probable values of 
/3^ can be determined by using the following criterion. 


= -^/?\A ^) 2 => minimum, 


where oo. is a statistical weighting function and is proportion to l/o- 2 , and 
a? is the statistical variance in the measurement in channel i . 

In formulating the linear least square analysis method we begin by assuming 
that the pulse height scale is invariant (i.e. , there is no gain shift or zero drift) 
and, therefore, the minimum can be found by taking a partial derivative with 
respect to some relative intensity . The derivative is then set equal to zero: 
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( 6 ) 


This can be written in matrix form as a no longer overdetermined set of equa- 
tions 

A cop - (kook) /3 = 0 , (6a) 

where the definition given in (4) still holds and 

co is a diagonal matrix of the weighting function a>.» and 
A is the transpose of the A matrix. Solving for the relative intensity vector 
/3, we find that 

fi = (kook) -1 A cop (6b) 

where (A«A) _1 is the inverse matrix of (kook) . Essentially it is the method of 
attaining the solution of equation (6b) which will concern us in the rest of the 
discussion. This equation has been applied to the solution of a number of spec- 
troscopic problems involving pulse height analyses 12, 13 > 14.15, 16, 17, 18 # There 
are a number of questions which must be answered in determining the algorithm 
(i.e., the method of solution) to be used for the given problem: 

1. Is it possible to invert the matrix (kook) ? 

2. Are the components of the A matrix linearly independent or do they 
correlate? 

3. How does one calculate this effect? 

4. How does resolution effect the correlation? 

5. What is the nature of the co matrix? 
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6. What is the effect of background subtraction? 

7. Is the system linear? 

8. If there are non-linearities in the system, how does one compensate for 
non-linearities in the application of this method? 

2. Library Functions: Let us consider the X-ray emission method, and 
determine the algorithm for obtaining /3 from equation (6b). We start by con- 
sidering the nature of the A matrix. The rows of the matrix can consist of 
either the monoelemental pulse height spectra (i.e., spectra characteristic of all 
the possible elements found in the irradiated sample) or the monoenergetic 
pulse height spectra (i.e., the pulse height spectra characteristic of all the 
possible monoenergetic x-rays striking the detector). Since we are interested 
in qualitative and quantitative, elemental composition, we will consider the prob- 
lem using the monoelemental pulse height spectra. 

The monoelemental pulse height spectra are determined by irradiating pure 
samples of the elements in the same geometric configuration as that used for 
the unknowns. Because these measurements have an inherent statistical error 
and because of the dynamics of the measurement system (possible gain drifts 
and zero drifts), these monoelemental pulse height spectra can be determined 
to only two or perhaps three significant figures with any great certainty. In 
the solution proposed in this paper, the statistical errors in the measurement 
of the monoelemental pulse height spectra are kept much smaller than the error 
in the measurement of the unknown pulse height spectra by extended and repeated 
measurements and therefore can be ignored. The difficulty comes then when an 
attempt is made to invert (Ao>A) . This inversion requires obtaining differences 
involving the third, fourth, and even fifth significant figure. As indicated above 
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the measurements of the monoelemental spectra are only good to two and at the 
most three significant figures. Thus, unless some additional constraints are 
introduced serious oscillations in the solution will result, and in fact one cannot 
successfully invert the matrix. Two necessary constraints that can be imposed 
are symmetry conditions of the (A«A) matrix and also a "non negativity option." 
The symmetry constraint follows from matrix algebra and is accomplished in 
the program which examines the symmetric off diagonal elements, i.e., the A 4j 
and A. . where i t j and sets them equal. The use of the non negativity option 
which proves very useful to the x-ray case will be further amplified. 

3. Non-negativity Constraint: It has been shown above that an exact solution 
cannot be obtained because of the statistical variation in the measurement. If 
such exact solutions were possible, the relative intensity of the various com- 
ponents would either be positive or zero. Now we are looking for the most 
probable solution, and in this case negative values can appear. When such nega- 
tive values are observed, there must be an associated error which is as great or 
greater than the absolute magnitude of the value obtained. Negative values with 
errors significantly smaller than the absolute value do however appear in the 
application of least square analysis to the reduction of pulse height spectra. It 
is this effect which produces oscillations in the solution described above. 

The oscillation in the solution can be described in the following way. A 
given relative intensity can appear as a negative value, and possibly more 
negative than it should be because of the errors in the inverse transformation 
just discussed. The overestimation of a negative component by the analytical 
method will tend to make some other component more positive to compensate. 
This in turn will tend to make a third component smaller than it should be in 
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order to compensate, and thus an oscillation in the solutions will be produced. By 
imposing a non-negativity constraint one essentially damps out these oscillations 
in the solution due to such factors as errors in the determination of the A matrix, 
and errors in improperly subtracting the background. What is meant by the non- 
negativity constraint then, is that when determining the relative intensities of the 
fi vector (6b), only positive values of monoelemental components are allowed, 
and negative values are set equal to zero. 

The use of physical constraints in the solution of linear least square analyses 
was first discussed by Beale 19 , and suggested for use in the analysis of pulse 
height spectra by Burrus 13 . A method using the non-negativity constraint for the 
analysis of gamma ray pulse height spectra 12 ’ 13 is now being applied to the prob- 
lem of x-ray non-dispersive analyses and has been discussed briefly in other 
papers 1,4 . 

Because of the necessarily empirical nature of the method of solution the 
absolute minimum obtained from equation (6b) is not necessarily the best solu- 
tion. One must find a minimum in the domain described by other physical con- 
straints. In addition to non- negativity, one should consider the following con- 
straints, for example: 

1. If the relative intensities of two species be known, then the ratio of in- 
tensities is invariant in the least square analysis. 

2. The sum of the relative intensities once determined for various elements 
always equals some constant. 

Later in the discussion concerning limits of resolution and correlation the 
latter condition will be invoked to help in the solution. Continuing with the 
discussion of a non-negativity constraint, we shall give only a brief and 
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non-rigorous physical argument to show how one obtains a solution (since a de- 
tailed description of this solution is given elsewhere 12 ’ 13 > 19 ).We first assume 
that we have a linearly independent set of monoelemental functions (function 
library A. . ) and we are required to find the relative intensities (/3^) of these 
library functions which combined will yield the best fit to the experimentally 
determined polyelemental pulse height spectrum (raw data spectrum p). The 
library must include all possible components. We now assume that all but two 
components from the library are zero. A least square fit (6b) to the total spec- 
trum p is used to determine the relative intensity of these two components. 
Since the number of components in the measured spectrum p is generally equal 
or larger than two, the estimates of the relative intensities of these two com- 
ponents will be greater than they should be. As components from the library 
are added, the relative intensities of the individual components will grow 
smaller. Thus if the least square fit gives a negative relative intensity for 
either of these two components, that component is eliminated from the library 
of elements and not used in the analyses (i.e., its relative intensity is set 
equal to zero). If both relative intensities are positive, both components are 
kept. Additional library components are added, one at a time, at least a square 
fit is made, and all negative components are again set equal to zero. The 
process is continued until all components have been tested. This then will yield 
the desired solution which is the least square solution for relative intensities in 
the positive domain. If the library elements are linearly independent then the 
solution obtained will be completely independent of the order in which the mono- 
elemental components are added. 
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The question may be asked concerning the necessity to perform the iterative 
process described above in order to find the solution to the minimum in the posi- 
tive domain. Why not, for example, find the absolute minimum with all the com- 
ponents at once and then adjust the negative values to zero. The answer is that 
because of the possibility of oscillation in the solution, it is difficult to tell 
whether the negative values result from the oscillations in the solution or be- 
cause of true statistical variation. The non-negativity, as we have indicated, 
tends to damp this oscillation out, and the algorithm described above finds the 
desired solution. 

One must be rather careful in using this non-negativity principal for strictly 
speaking, negative solutions should be allowed in least square analysis. What 
the non-negativity constraint does is to increase the mathematical error because 
one does not find the absolute minimum as required by the least square method. 
Furthermore, by rejecting certain components, estimates of possible minimum 
detectable limits for these components are not determined. In order to obtain 
the minimum detectable limits in the program prepared for this analysis, both 
the absolute minimum and the minimum in the positive domain can be found, and 
compared. 

4. Error Calculation : We now consider the problems of calculating the error 
in the relative intensity /3^ due to the statistical variation in p . , the counts in 
channel i , contributed by the components of the polyelemental pulse height 
spectrum. Let us consider equation (6b) , which shows that can be written as 
a linear combination of the p i ’s: 
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We define C = (Ao>A), where (A^A) is a symmetric matrix whose elements 


of C are given by 


~vy 




Remembering that a), is the statistical weight and a>. = 1 /a\ where cr\ is the 
variance of the measurement of p { in channel i, and C _1 is the inverse of 


matrix C, e.g.. 


CC" 1 = I, 


( 9 ) 


and where I is the identity matrix (i.e. , all off diagonal elements are zero and 
the diagonal elements are one) or 


U = 1 if * = * 
Kk = 0 if v * x 


( 10 ) 


The elements of the can be given 

U = J2 c -y c *~' (11) 

y 

From equation (6c), it is seen that /3^ is a linear homogeneous function of the 
counts p . under the assumption that there is no significant error in the A^ 
compared to the measurement of p . . Thus the mean square deviation a 2 (f3^) 
corresponding to the variation in p can be written as a linear sum of the vari- 
ances on the Pj’s. Thus from equation (6c) we get the variance for /3 to be 

°- 2 w = Z] H H c ^~ 1 1 \ \ ^ a2 (Pi ) ’ (l2> 

ivy 
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and from (8), 


= 



CO. 

l 


or 



v K 


L 


CO. 


A. 


Then from (7), 



'■W = C A _1 C *- 1 c vr 

v y 


or 


v y 


and from equation (11) 


: (/\> = £ 


C v \ _1 x v\- 


Finally from (10) 


03 (^\) = c ^- 1 , 


(13) 


that is, a 2 (/S^_) can be found from the diagonal elements of the C~ 1 matrix. 
Equation (13) is true if yi 2 is equal to one. Xi 2 is defined as follows: 


xf = 


cp, - z 

i A. 


Pk A iX ) 2 


n - m 


( 12 ) 


where n-m are the number of degrees of freedom, 
n is the number of channels, and 
m is the number of library components. 
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If X i t 1, then 


0x> = x? c KK 


-X 


(13a) 


xf and its utilization will be discussed in slightly more detail below. 

5. Correlation : Let us now consider what we mean by linearly independent 
library functions. It was shown that C^ - 1 is the variance on fix. It can be 
further shown that C^" 1 is the covariance of A.'th and y'th component 20 ’ 21 . 

A measure of interferences F^ y between the A/th and y'th component is given 

by 








100 % 


(14) 


Equation (14) is a measure of how different, or how well one component or library 
spectrum can be resolved from another. It is also a measure of whether or not 
the components in the library of spectra can be considered to be linearly inde- 
pendent. 

As an illustration of how equation (14) is used, let us assume that the library 
function or monoelemental functions can be described by Gaussians: 


= exp 


_ (i - P^) 2 


2a 


k 


(15) 


where is the normalized counts in channel i due to the A-'th component, 
a^ is a constant and is a measure of the width of the Gaussian for energy k, 
and was defined in equation (1). 

Now let us calculate the percent inference between two monoenergetic pulse 
height spectra with Gaussian form using equation (14) and a shape given by (15). 
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If we assume no statistical error, equation (14) can be written 


E ^ 




(14a) 


Using equations of the form given in (15) for A iy and A lk and replacing the sum- 
mation by an integral from - infinity to + infinity (14a) becomes 


yK 


/* +00 

(i -P.J 2 

exp - exp - 

^ - Py) 2 

*2-00 

2a 2 

H 


di 


*+co 


(di) exp - 


(i -P 7 ) 2 


n2 




.+<» 


(i - p k y 

exp (di) 

2a l 


-1 2 


(14b) 


2a \ a y 
a £ + a y 


exp - 


Py-Py 

‘ a £ + a 2 


The constant a y can be written in terms of the resolution R y as 




2 / 2 Ln 2 


(16) 


Remember that P y stands for the pulse height portion of some energy E y and is 
proportional to E y . If we now assume that equation (2) is applicable and let 
n = . 5, we get 


R y 


Py 

P r 


(17) 


For this special case we can substitute equation (17) into equation (14b) 


■yX. 


2 /PyP 

Py + P, 


y Py (P\ - Py ) 2 


exp - 


(^- p rW 


(14c) 
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Figure 7 is a plot of Equation (14c) for various detector resolutions as a function 
of the percent separation (P^- P^)/P^for decreasing P y , where P y and are 
the axis for two adjacent energies. If P y = P\, the Gaussians are identical and 
the percent interference is 100%. As (P y - P\ ) increases the percent interfer- 
ence decreases, and this is then a measure of how well the two Gaussians can 
be resolved. For example, consider the following case. Given 

1. Two x-ray energies separated by 10%, 

2. The detector resolution for the higher energy is 20%, 

3. The relationship given in equation (17) holds, and 

4. The photopeak is described by a Gaussian which fully described the pulse 
height spectrum of the monoelemental component (e.g., no escape peak or 
continuum), there will be a 50% interference between these energies. 

The sum fi y + /3^ can be determined and will be constant, but both fi y and 
/. 3 ^ can be varied by 50% keeping their sum constant and the fit based on the least 
square criteria will be just as good. In general given a library of functions A, 
the percent interference between various components can be determined by 
forming (AA) -1 and calculating F ^ from (14). This point will be described 
further in Section III, where a solution of an actual experimental problem is 
presented. 

If two functions do interfere strongly, it is sometimes possible to impose a 
physical constraint, which can eliminate the interference. For instance, in the 
problem of gamma ray spectroscopy with a mixture of isotopes with varying 
half lives, the pulse height spectra can be followed until one of the species has 
decayed out. The second specie is then followed and the intensity determined. 
Knowing the half life of the second species, the intensity of the second component 
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can be determined and extrapolated back to the earlier measurements when the 
correlation was large. The second component can then be stripped out and the 
analysis repeated for the first element with the correlation due to the second 
element eliminated. In x-ray fluorescence, on the other hand, it is possible to 
use filters, and determine either the intensity of one energy or determine the 
ratio of one energy with respect to another. Using the intensity of one component 
or the ratio of one component to another as a physical constraint, the interference 
can be eliminated. 

The percent interference described above has only considered the problem 
of resolution, but there are two other effects which will produce strong inter- 
ference. In one case the statistical error due to counting and background sub- 
traction may be so high that the percent interference will be large. This will be 
related to the statistical weight a> in the (AcoA)" 1 equation. Secondly, interference 
can occur between two elements if the library is incomplete and an element is 
missing. Then the two components on either side of the missing element will be 
under or overestimated in order to make up for the missing element and the two 
adajcent elements will correlate. When one calculates the difference between the 
measured spectrum and the calculated spectrum using least squares, a negative 
or positive peak will appear at the position of the missing element. By adding in 
the missing component, the interference and x \ can be greatly reduced. 

6. Chi Square: In the method developed for this problem ax? test is per- 
formed (see Equation 12) for goodness of fit. The closer X \ is to unity the 
better the fit. 

7. Background Correction: At first thought, this problem should be rather 
simple to handle. For example, a background is measured for either the same 
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time or a time longer than the measurement of the unknown system, and then the 
background is subtracted with adjustment being made for the counting times. In 
x-ray fluorescence a major part of the background is due to coherent scattering 
from the sample of the incident radiation. Thus, the fraction of background to be 
subtracted will be affected by the nature of the sample (average atomic number, 
density, etc.). However, as an alternative approach the following method was 
used. Measurements of the raw data spectrum p are made to extend to energies 
(or pulse height channels) higher than the highest energy expected in the pulse 
height spectrum due to the sample. A background spectrum using a strong 
scatterer such as boric acid or plastic in the sample position is taken. Charac- 
teristic x-rays from these substances are so soft as to remain undetected by our 
present instrumentation. As was pointed out in the section on the shape of the 
pulse height spectra, the shape of the background spectrum does not seem to be 
strongly affected by the scatterer although the intensity is greatly affected. The 
background is then included as a library component and subjected to the least 
square treatment, already described, which yields a measure of the background 
intensity included in the total spectrum. In practice this latter approach has 
proven feasible because in the measurement of the raw data spectrum p there 
is a higher energy portion of the background spectrum which is free of lines 
generated in the specimen. If this were not the case, the background spectrum 
would strongly correlate with each of the monoelemental components, and a 
unique solution would not be obtained. In this event one goes back to the sub- 
traction technique. This second approach has also been used, but since an 
iterative process is necessary, it is more time consuming. This follows because 
the fraction of background to be subtracted is not known and can only be arrived 
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at by trial and error, e.g., various fractions are subtracted, subjected to least 
square analysis and a minimum y? sought. In this study, the results presented 
in the paper used the method of including the background as a library component 
although the computer program developed has the subtract option included. In- 
cluded as part of this second option is an error calculation due to this subtract 
mode and a proper evaluation of the weighting function co used in the least square 
analysis. 

8. Gain Shift: Up to this point it has been assumed that the pulse height 
scale for the library spectra is the same as that for the measured raw data 
spectrum p. Because of the possibility of gain shifts, and zero drift, the pulse 
height scale can be compressed or expanded and linearly shifted. Two approaches 
can be used to compensate for this effect. Non-linear least square analysis 
methods may be used. That is, in equation (5), the derivative with respect to 
yQ and A. y can be obtained. This technique is rather complex, but a number of 
methods have been developed 22, 23 . Because of these complexities we have used 
the linear method, and make corrections to the raw data vector p. before per- 
forming the analysis. Other similar techniques have been developed 13,15,18,24 . 

Let us first consider the problem of gain shift. It must be remembered that 
the measured pulse height spectrum is a histogram, that is, it is a measure of 
all pulses in some increment Ai about i as a function of i . The total sum £ p i 

i 

is equal to the total number of interactions between the indicent x-ray flux and 
detector. This number must remain a constant for a given measurement and is 
independent of a compression, extension, or linear translation along the pulse 
height axis. If we assume that there is a gain shift G on the pulse height axis 
between library spectra A. . and the raw data spectrum p , the pulse scale i 
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will have to be multiplied by G and the count rate scale p. divided by G to com- 
pensate for the shift. This keeps the area a constant. 

It is important to point out that the library spectra histograms are only 
included for integer values of pulse height which are separated by Ai = 1. The 
gain is shifted by some value G which will produce values at fractional pulse 
heights and cause Ai to become either greater or less than unity. For example 
consider the case for G = .98 or G = 2.0. Channels ten and eleven will become 
channels 9.8 and 10.78 for G = .98 and channels 20 and 22 for G = 2.0. For the 
first case G - .98, we must find the value at channel 10, 11, 12, etc. and for the 
second G = 2.0 we would have to find the value at channel 20, 21, 22, etc. The 
value at channel 10 (G = .98) or at channel 21 (G = 2.0) is found in our method by 
linear extrapolation which is an integral part of our program. Only two points 
are used at a time in this extrapolation. Higher order methods can also be used, 
but have not been found necessary in this application. 

Two choices are possible: shifting the library spectra to correspond to 

the raw spectrum or visa versa. The second alternative has been found more 
economical in terms of time. The results obtained using either alternative 
are the same. 

In practice either a given gain shift or a range in gain shift (G min and G max ) 
is assumed. A series of least square fits between the assumed G . and G 

min max 

are made, y 2 calculated and the minimum y 2 found for this range. Figure 8 
shows a measured spectrum and the gain shifted spectrum. The spectrum has 
been shifted 30%. 

A similar method compensating for zero drift is now being developed for 
inclusion in the program. Here a minimum e and maximum e is 

min max 
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introduced. The pulse height scale is shifted by a value e, the intercepts at 
integer pulse height values are then determined, and again least square fits are 
made until a minimum y? is determined in the range e . to e 

9. Computer Program: The computer program flow diagram which has been 
developed to perform the above described calculations is shown in Figure 9. The 
zero shift part of the program is not included in the flow diagram. There are 
four flows available in this program. Flow 1 allows the stripping of one or more 
components from the spectrum before analysis. Flow 2 allows background sub- 
traction, gain shifting of the spectrum, and contains the non-negative constraint 
as an option. Flow 3 is just a flow for solving equation (5) without gain shift and 
without subtraction of background, and without the non-negativity constraint. 

Flow 4 is similar to Flow 2, but does not allow for background subtraction. 
Calculation for errors introduced due to stripping of components and background 
subtraction are included in the appropriate flow. A number of options are al- 
lowed for calculation of oo • w, the statistical weight matrix, can be either (1) set 
equal to the unity matrix, (2) set equal to 1/p, (3) calculated internally in Flows 
1, 2 and 4, or (4) read in separately according to the analysts desire. The output 
from the computer program yields the following information: (1) the flow used, 
(2) the gain for minimum y? if the option is chosen, (3) whether or not the non- 
negativity constraint was used, (4) the relative intensities for /3 , (5) the mono- 
elemental library functions used, (6) the elements rejected if non-negativity is 
used, (7) the standard deviation corresponding to each of the relative intensities 
calculated, (8) the percent interference between library elements, (9) the actual 

c 

spectrum analyzed, (10) the synthesized spectrum obtained using the relative 
intensities determined by least square, (11) the difference between the actual and 
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Figure 9 
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synthesized spectrum, (12) and the value. A more detailed description of the 
program can be found in references 10 ’ n • 22 . 

ID. APPLICATION OF LEAST SQUARE TECHNIQUE 
TO AN EXPERIMENTAL PROBLEM 
A. General Experimental Results 

The procedures which have been described in some detail above have been 
applied to the problem of the semi-quantitative determination of the composition 
of a variety of rock types ranging from ultra-basic rocks such as dunite and 
peridotite through acidic rocks such as granite. For this purpose a suite of six 
rocks, carefully prepared by the U. S. Geological Survey and presently being 
circulated for comparative analysis were used. All rock samples were finely 
divided and homogenized. 

Library spectra were obtained using oxides and carbonates of the various 
elements to be determined such as Mg, Al, Si, K, Ca, Ti, Mn and Fe. 

For the measurements reported here, a sealed A-CH 4 proportional counter 
was used which gave a resolution of 17.5 percent of the Mn Ka line from Fe ss . 

Background was determined using a boric acid briquet as a scatterer. 
Figures 10-15 show the observed spectra for dunite, peridotite, basalt, andesite, 
granodiorite and granite. These experimentally measured data are shown as 
points superimposed on the solid line curves, synthesized by the computer 
using the least square technique proposed and the library elements. 

In every instance the spectra have been gain shifted so that they are on a 
common pulse height energy scale. On this basis one can state that up to ap- 
proximately channel 90, covering the elemental range from Mg to Fe, the agree- 
ments are excellent and in fact within the statistical error. Beyond channel 90 
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where one observes a number of descrete lines in the background spectrum, 
there is a slight divergence, very likely due to the difference in scattering be- 
tween the samples measured and the boric acid pellet. These divergences shown 
in Figures 10-15 would actually appear smaller if the statistical spread for the 
data points were plotted (see Figure 6). Nevertheless these differences are real 
and do increase the statistical error in the determination of the relative intensi- 
ties. These increased errors are reflected in an increase of x? . The increase 
on the variance a 2 (j3^) due to increase in y 2 is then given by equation (13a) 
above. 

In order to produce the synthesized curves shown in Figures 10-15, the 
relative intensities of each of the library components has been determined by 
using the least square criteria and the non-negativity constraint. These relative 
intensities are shown in Table I together with the statistical error. Blank por- 
tions in the table indicate components showing negative or zero intensity values. 
These are automatically rejected in the analysis and not included as part of the 
library. This is an indication that a particular component is either absent or its 
intensity contribution is indistinguishable statistically from the noise. The 
actual composition for each sample is shown in the corresponding figure (10-15). 
The percentage interference between elements (i.e., those greater than 10%) are 
given in Table II. It can be seen that the elements fall into three groups where 
interference is strong. These are: (1) Si, Mg, Al, (2) Ca, K, and (3) Fe, Mn. 

In these groups, unless an additional constraint is used such as that described 
below, it is very difficult to determine the concentrations of the individual 
components even though the sums of the relative intensities are correct. 
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Table I 

Relative Intensities Obtained from the Least Squares Analysis 
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Table H 


Correlation Greater Than 10% between 
Library Components 



Si 

A1 

Fe 

Mg 

Ca 

K 

Ti 

Mn 
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B. Resolving Correlations 

We now shall consider the case of the Si-Al-Mg correlation in detail. As 
indicated above, another physical constraint is required in order to resolve the 
interferences. Since aluminum has a strong absorption edge lying between the 
silicon and the aluminum K lines, a 0.5 mil aluminum foil was used in front of 
the detector window to preferentially absorb the silicon K radiation. Figures 
16 and 17 show for example the pulse height spectra of the dunite and granite 
taken under these conditions. In comparison to the spectra for these rocks 
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Table HI 


Results Obtained for Plan View Meteorite Using the Least Squares Analysis 


Element 

Relative Intensities from 
Least Squares Analysis 

% Composition 
from Fig. 17 

Reported 

Chemical Analysis* 

Fe 

.36 x 10 s ± .02 x 10 s 

14.1 

± 2.8 

17.-24. 

Si 

.31 x 105 ± .01 x 105 

17.5 

± .7 

17.25 

A1 

.25 x 10 4 ± .09 x 10 4 

3.2 

± 1.5 

1.08 

Mg 

.54 x 10 4 ± .05 x 10 4 

18. 

± 2.5 

13.71 

Ca 

.24 x 10 4 ± .03 x 10 4 

1.08 ± .21 

1.19 

K 

.95 x 103 ± .30 x 103 

.1 

± .03 

.066 


♦Quarterly Rept of the U.S. Geological Survey, April-June 1965. 

previously shown in Figures 10 and 15 the dunite now shows a distinct Mg peak 
and the A1 line is predominant in the granite. 

The least square analysis technique using a non-negativity constraint and an 
aluminum foil filter was applied to the same suite of six rocks shown in Figures 
10-15. Because of the strong silicon absorption, the silicon component was 
eliminated from the analysis. The relative intensities of the A1 and Mg obtained 
were corrected for absorption due to the A1 filter. While a strong correlation is 
possible between these elements, this was not a significant factor in this instance 
because where the Mg concentration was high, the A1 concentration was very 
low and vice-versa. 
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If A1 and Mg were both present in large amounts or if it is required that the 
A1 and Mg be known with great accuracy, then an magnesium filter could be used 
to selectively absorp the A1 radiation, and resolve this correlation. 

Since the A1 and Mg relative intensities had been determined, these com- 
ponents were then stripped out of the original spectrum and the Si relative in- 
tensity determined using Flow 1 of the computer program. These results are 
listed in Table I as corrected values for Si, A1 and Mg and can be compared to 
the uncorrected values within the same block. 

Figure 18 is a plot on a log- log scale of the relative intensities versus chemi- 
cal composition for Si, Al,Mg,K,Ca and Fe. The solid lines are least square fits 
to the data points obtained from the computer analysis. There is a good approxi- 
mation to a linear relationship between relative intensity and percentage composi- 
tion, consistent with the small concentration range. No line was drawn for the Mg 
because the lower Mg concentrations are associated with high A1 concentrations in 
these samples and the strong correlation between these two elements makes it 
difficult to determine the Mg component without very large errors. It is significant 
that one can easily distinguish between a 10% relative concentrational variation in 
the 30% range. The error bars in these plots are strictly fluctuations due to counting 
statistics and background correction. These can be partially reduced by increased 
counting times and background reduction. 

C. Classification of Rock Types 

Although the curves in Figure 18 show a reasonable linear approximation 
between relative intensity and concentration one does not intend to infer that this 
method can yield precise chemical analysis. No attempt has been made at care- 
ful sample preparation or to correct for matrix effects such as absorption or 
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enhancement. In fact, in terms of a lunar surface mission, a serious attempt 
has been made to determine how much chemical information can be obtained 
without careful control above factor. Using these ground rules the data will be 
used in an attempt at broad classification of rock types. The approach is 
similar to those in the studies performed using neutron methods described in 
refernces 25-29. 

A preliminary classification of the rocks according to type is made on the 
basis of the presence or absence of certain chemical elements simply determined 
from the computer output, for example the presence or absence of large amounts 
of Mg, Al, Ca and K. For example the ultra-basic rocks are rich in Mg with 
negligible amounts of Al, Ca and K (see Table I and Figures 10-15). 

Following this, an approach similar to that described for the neutron-gamma 
techniques (references 25-29) is used. Rather than attempting to carefully control 
the sample geometry, ratios of various elements have been calculated as shown 
in Figures 19-25. One observes again a good approximation to a linear relation- 
ship between chemically determined and calculated ratios. As pointed out by 
Waggoner (reference 29) it is obvious that no one element or ratio of elements can 
be used to classify a rock type. However, if one combines a number of these 
factors, then it is possible to categorize a rock as lying somewhere on the scale 
between the ultra-basic and acidic rocks. This can, of course, only be done to 
the extent that a chemical analysis is meaningful in this context. 

D. Quantitative Analysis 

The performance of quantitative analysis normally requires extreme care 
in sample preparation, the use of appropriate standards and adjustment for ma- 
trix effects. Such requirements are sometimes difficult to fulfill even in a 
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laboratory environment. One cannot even hope to control these parameters in a 
remote exploration program. Keeping these constraints in mind, an attempt was 
made, using these rough procedures to perform an analysis on a meteorite to 
obtain some idea of the possibility for performing a quantitative analysis. 

Figure 26 shows the measured pulse height spectrum obtained for the Plainview 
Meteorite, a bronzite-chondrite. Table III summarizes the results obtained 
using the curves of Figure 17 as calibration curves. The agreements can be 
considered as quite satisfactory for a first cut chemical analysis. These results 
combined with the data used to prepare the calibration curves show that one can 
expect to obtain useful chemical values, satisfying the goals for a lunar geo- 
chemical exploration device. 
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LIST OF ILLUSTRATIONS 


1. a source x-ray analyzer 

2. a emitter source holder 

3. Background pulse height spectrum: cm 242 source; boric acid scatterer; a 
90% A, 10% CH 4 , .001" Be window proportional counter. 

4. Counter tube efficiency: 90% A, 10% CH 4 filling; .001 in Be window; 2 cm 
absorption path. 

5. Discrete energy spectrum and corresponding pulse spectrum of the char- 
acteristic Fe x-ray lines. Pulse height spectrum measured with a sealed 
90% A, 10% CH 4 proportional counter. 

6. Pulse height spectrum of basalt measured with a proportional counter filled 
with 90A and 10% CH 4 - The following can be seen: (1) The monoelemental 
components; (2) The envelope due to the sum of the monoelemental compo- 
nents; and (3) The measured pulse height spectrum with statistical variations 
included. 

7. Percentage of interference for various resolutions. 

8. The effect of gain shift on silty sand from Hoppe Butte, Arizona. 

9. Flow diagrams for least squares analysis. 

10. Pulse height spectrum of dunite. 
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11. Pulse height spectrum of peridotite. 

12. Pulse height spectrum of Basalt. 

13. Pulse height spectrum of andesite. 

14. Pulse height spectrum of grandiorite. 

15. Pulse height spectrum of granite. 

16. Dunite spectrum with l/4 mil aluminum absorber. 

17. Granite spectrum with 1/2 mil aluminum absorber. 

18. Measured relative intensity as a function of chemical composition. 

19. Measured Fe/Si intensity ratios versus chemical composition. 

20. Measured Al/Si intensity ratios versus chemcial composition. 

21. Measured Mg/ Si intensity ratios versus chemical composition. 

22. Measured Cu/Si intensity ratios versus chemical composition. 

23. Measured K/Si intensity ratios versus chemical composition. 

24. Measured (Cu + K)/Si intensity ratios versus chemical composition. 

25. Measured K/Cu intensity versus chemical composition. 

26. Spectrum of Planvie Meteorite (bronze). 
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